/*******************************************************************************
Project : Bank Expansion and Moneylender Interest Rates - RDD Evidence from India
Authors : Briones L., Chiseliov V., Gushchin E., Guttensen J.
File    : All Figures.do
Purpose : Master script to generate all figures in the replication package.
*******************************************************************************/

******************************************************
* SETUP AND DIRECTORY
******************************************************
cd "`c(pwd)'\Data"

local figure_01 1   // Done		// Map
local figure_02 1   // Done 	// hist_cutoff
local figure_03 1   // Done		// int hist
local figure_04 1   // Done		// fuzzy RD
local figure_05 1   // Done		// smoothness
local figure_08 1   // Done		// Heterogeneity

if `figure_01' {

* Install required packages
ssc install spmap, replace
ssc install shp2dta, replace

clear all

******************************************************
* SHAPEFILE DOWNLOAD AND CONVERSION
******************************************************
local shp_url "https://github.com/datameet/maps/raw/master/Districts/Census_2001/2001_Dist.shp"
local dbf_url "https://github.com/datameet/maps/raw/master/Districts/Census_2001/2001_Dist.dbf"
local shx_url "https://github.com/datameet/maps/raw/master/Districts/Census_2001/2001_Dist.shx"

local shp_file "Shapefiles/2001_Dist.shp"
local dbf_file "Shapefiles/2001_Dist.dbf"
local shx_file "Shapefiles/2001_Dist.shx"

copy "`shp_url'" "`shp_file'", replace
copy "`dbf_url'" "`dbf_file'", replace
copy "`shx_url'" "`shx_file'", replace

shp2dta using "`shp_file'", ///
database("Shapefiles/india_districts.dta") ///
coordinates("Shapefiles/india_coords.dta") ///
genid(A_MapID) replace 

use "Shapefiles/india_districts.dta", clear

******************************************************
* CLEAN AND STANDARDIZE NAMES
******************************************************
gen State_UT = ST_NM
gen District = DISTRICT 

replace State_UT = "Andaman & Nicobar Islands" if ST_NM == "Andaman & Nicobar Island"
replace State_UT = "Arunachal Pradesh" if ST_NM == "Arunanchal Pradesh"
replace State_UT = "Dadra & Nagar Haveli" if ST_NM == "Dadara & Nagar Havelli"
replace State_UT = "Jammu & Kashmir" if ST_NM == "Jammu and Kashmir"
replace State_UT = "Manipur +" if ST_NM == "Manipur"
replace State_UT = "Pondicherry" if ST_NM == "Puducherry"
replace State_UT = "Uttaranchal" if ST_NM == "Uttarakhand"
replace State_UT = "Delhi" if ST_NM == "Delhi  & NCR"

replace District = "New Delhi" if ST_NM == "Delhi  & NCR"
replace District = "Greater Mumbai" if DISTRICT == "Mumbai" | DISTRICT == "Mumbai (Suburban)"
replace District = "Perambalur" if DISTRICT == "Ariyalur"
replace District = "Hazaribag" if DISTRICT == "Hazaribagh"
replace District = "Senapati" if DISTRICT == "Senapati (Excl. 3 sub-divisions)"
replace District = "Sant Ravidas Nagar (Bhadohi)" if DISTRICT == "Sant Ravidas Nagar Bhadohi"
drop if District == "Data Not Available"

******************************************************
* LINK TO CENSUS_UB DATA
******************************************************
reclink State_UT District ///
using "MOF_Data\Census 2001\Census_UB.dta", ///
idmaster(A_MapID) idusing(ID_Census) ///
gen(match_score)

save "Shapefiles/india_final.dta", replace

******************************************************
* CREATE UNDERBANKED MAP
******************************************************
use "Shapefiles/india_final.dta", clear

label define underbanked_labels_new 0 "Not Treated" 1 "Treated"
label values Underbanked underbanked_labels_new

spmap Underbanked using "Shapefiles/india_coords.dta", ///
id(A_MapID) ///
fcolor("117 101 143" "235 179 076") ///
clmethod(unique) ///
legstyle(3) ///
leglabel("Not Treated" "Treated") ///
ndocolor(none) ///
ndlabel("No data") ///
legend(region(lstyle(none)) position(1) size(small)) 

graph export "output/figures/heatmap_underbanked.pdf", replace
graph export "output/figures/heatmap_underbanked.eps", replace
}


if `figure_02' {
***** histogram ******
use "MOF_Data\RBI_2005_Q1\MOF_Census_2005_Q1.dta", clear
drop if B_Ratio_Dist_Ex == . // which ones are deleted?

hist ratio_banks if ratio_banks<20000 & ratio_banks>10000, freq width(100) fcolor("80 106 136%75") ///
		xline(15070, lpattern(dash) lcolor("217 85 130")) lcolor("80 106 136%75") ///
		xline(11183, lpattern(dash_dot) lcolor(black)) xline(18957, lpattern(dash_dot) lcolor(black)) ///
		xtitle("District level population-to-branch ratio", size(large)) ///
		ytitle(, size(large)) ylabel(, labsize(large)) xlabel(,labsize(large)) ///
		plotregion(fcolor(white)) graphregion(margin(l=10 r=10) fcolor(white) lwidth(0)) ///
		
graph export "output/figures/histogram.pdf", replace
graph export "output/figures/histogram.eps", replace
}


if `figure_03' {

*-------------------------------*
* Sample I
*-------------------------------*
use AllIndiaI, clear

* Remove outliers and missing values
keep if abs(stdlog_i) < 4
drop if missing(log_i, B_Ratio_Dist_Ex)

* Histogram
hist log_i,  width(0.28) color("37 150 190") graphregion(color(white)) plotregion(color(white))
	
* Export
graph export "output/figures/irate_wavei.png", replace
graph export "output/figures/irate_wavei.eps", replace

*-------------------------------*
* Sample II
*-------------------------------*
use AllIndiaII, clear

* Remove outliers and missing values
keep if abs(stdlog_i) < 4
drop if missing(log_i, B_Ratio_Dist_Ex)

* Histogram
hist log_i, width(0.28) color("37 150 190") graphregion(color(white)) plotregion(color(white))
	
* Export
graph export "output/figures/irate_waveii.png", replace
graph export "output/figures/irate_waveii.eps", replace

}


if `figure_04' {
*** threshold
use "MOF_Data\RBI_2005_Q1\MOF_Census_2005_Q1.dta", clear

drop if B_Ratio_Dist_Ex == . // which ones are deleted?

scatter Underbanked B_Ratio_Dist_Ex if B_Ratio_Dist_Ex<20000, xlabel(,labsize(large)) ///
			plotregion(fcolor(white)) graphregion(margin(l=10 r=10) fcolor(white) lwidth(0)) ///
			xline(0, lpattern(dash) lcolor("217 85 130")) ytitle("Treatment status", size(large)) ///
			xtitle("District level population-to-branch ratio (normalized)", size(large)) ///
			ylabel(0 1, labsize(large)) msize(small) mcolor("235 179 076%20") ///
		
graph export "output/figures/fuzziness.pdf", replace
graph export "output/figures/fuzziness.eps", replace
}


if `figure_05' {
/****************************************************/
/* SMOOTHNESS CHECKS: CREDIT MARKET CHARACTERISTICS */
/****************************************************/

/* SMOOTHNESS CHECKS: INTEREST RATE FROM MONEYLENDERS */

use AllIndiaI, clear

global h = 7126

/* Keep only moneylender loans */
keep if inlist(B15_2_q6,"10","11") /*Moneylender*/ 

/* Remove outliers */
keep if abs(stdlog_i)<4  /* remove observations beyond 3 standard deviations (2.5%)*/
drop if log_i==.
drop if B_Ratio_Dist_Ex==.

destring B15_2_q3 B15_2_q6, replace

duplicates drop HHID, force

rdplot log_i B_Ratio_Dist_Ex if inlist(B15_2_q3, 2000,2001,2002,2003) & inrange(B_Ratio_Dist_Ex,-$h,$h), p(1) weights(Weight) binselect(esmv) fuzzy(treatment) graph_options(graphregion(color(white))) genvars

/* Create custom plot */
set scheme s2color
tw (lfitci log_i  B_Ratio_Dist_Ex if B_Ratio_Dist_Ex<0 & B_Ratio_Dist_Ex>-$h) ///
   (lfitci log_i  B_Ratio_Dist_Ex if B_Ratio_Dist_Ex>0 & B_Ratio_Dist_Ex<$h), ///
   xline(0, lpattern(dash) lcolor("217 85 130")) ///
   legend(off) ///
   xtitle("District level population-to-branch ratio (normalized)", size(large)) ///
   xlabel(, labsize(large)) /// Increase font size of X-axis tick labels
   ylabel(, labsize(large)) /// Increase font size of Y-axis tick labels
   plotregion(fcolor(white)) ///
   graphregion(fcolor(white) lwidth(0)) || (scatter rdplot_mean_y rdplot_mean_x, msize(small) mcolor("45 152 165"))
   
graph export "output/figures/smoothness_interestrate1.pdf", replace
graph export "output/figures/smoothness_interestrate1.eps", replace


/* SMOOTHNESS CHECKS: LOAN AMOUNT FROM ALL SOURCES */

use AllIndiaI, clear

global h = 7126

destring B15_2_q3 B15_2_q6, replace

rename B15_2_q3 b14_q3

/* Merge CPI data for inflation adjustment */
merge m:1 b14_q3 using CPI_index, keepusing(CPI_index)
drop if _merge != 3
gen loan_amount = log(B15_2_q5 /(CPI_index * 10))

rdplot loan_amount B_Ratio_Dist_Ex if inlist(b14_q3, 2000,2001,2002,2003) & inrange(B_Ratio_Dist_Ex,-$h,$h), p(1) weights(Weight) binselect(esmv) fuzzy(treatment) graph_options(graphregion(color(white))) genvars

/* Create custom plot */
set scheme s2color
tw (lfitci loan_amount  B_Ratio_Dist_Ex if B_Ratio_Dist_Ex<0 & B_Ratio_Dist_Ex>-$h) ///
   (lfitci loan_amount  B_Ratio_Dist_Ex if B_Ratio_Dist_Ex>0 & B_Ratio_Dist_Ex<$h), ///
   xline(0, lpattern(dash) lcolor("217 85 130")) ///
   legend(off) ///
   xtitle("District level population-to-branch ratio (normalized)", size(large)) ///
   xlabel(, labsize(large)) /// Increase font size of X-axis tick labels
   ylabel(, labsize(large)) /// Increase font size of Y-axis tick labels
   plotregion(fcolor(white)) ///
   graphregion(fcolor(white) lwidth(0)) || (scatter rdplot_mean_y rdplot_mean_x, msize(small) mcolor("45 152 165"))
   
graph export "output/figures/smoothness_amountborrowed1.pdf", replace
graph export "output/figures/smoothness_amountborrowed1.eps", replace

/* SMOOTHNESS CHECKS: PROB OF TAKING BANK LOAN */

use AllIndiaI, clear

bysort HHID: gen N = _N

rdplot N B_Ratio_Dist_Ex if inrange(B_Ratio_Dist_Ex,-$h,$h), p(1) weights(Weight) binselect(esmv) fuzzy(treatment) graph_options(graphregion(color(white))) genvars

/* Create custom plot */
set scheme s2color
tw (lfitci N  B_Ratio_Dist_Ex if B_Ratio_Dist_Ex<0 & B_Ratio_Dist_Ex>-$h) ///
   (lfitci N  B_Ratio_Dist_Ex if B_Ratio_Dist_Ex>0 & B_Ratio_Dist_Ex<$h), ///
   xline(0, lpattern(dash) lcolor("217 85 130")) ///
   legend(off) ///
   xtitle("District level population-to-branch ratio (normalized)", size(large)) ///
   xlabel(, labsize(large)) /// Increase font size of X-axis tick labels
   ylabel(, labsize(large)) /// Increase font size of Y-axis tick labels
   plotregion(fcolor(white)) ///
   graphregion(fcolor(white) lwidth(0)) || (scatter rdplot_mean_y rdplot_mean_x, msize(small) mcolor("45 152 165"))
   
graph export "output/figures/smoothness_numberofloans1.pdf", replace
graph export "output/figures/smoothness_numberofloans1.eps", replace

/* SMOOTHNESS CHECKS: PROB OF TAKING BANK LOAN */

use AllIndiaI, clear

global h = 7126

destring B15_2_q3 B15_2_q6, replace

/* Create credit agency variables */
gen banks = (inlist(B15_2_q6, 2, 3))
replace banks = . if B15_2_q6 == .

duplicates drop HHID, force

rdplot banks B_Ratio_Dist_Ex if inlist(B15_2_q3, 2000,2001,2002,2003) & inrange(B_Ratio_Dist_Ex,-$h,$h), p(1) weights(Weight) binselect(esmv) fuzzy(treatment) graph_options(graphregion(color(white))) genvars

/* Create custom plot */
set scheme s2color
tw (lfitci banks  B_Ratio_Dist_Ex if B_Ratio_Dist_Ex<0 & B_Ratio_Dist_Ex>-$h) ///
   (lfitci banks  B_Ratio_Dist_Ex if B_Ratio_Dist_Ex>0 & B_Ratio_Dist_Ex<$h), ///
   xline(0, lpattern(dash) lcolor("217 85 130")) ///
   legend(off) ///
   xtitle("District level population-to-branch ratio (normalized)", size(large)) ///
   xlabel(, labsize(large)) /// Increase font size of X-axis tick labels
   ylabel(, labsize(large)) /// Increase font size of Y-axis tick labels
   plotregion(fcolor(white)) ///
   graphregion(fcolor(white) lwidth(0)) || (scatter rdplot_mean_y rdplot_mean_x, msize(small) mcolor("45 152 165"))
   
graph export "output/figures/smoothness_bankloan1.pdf", replace
graph export "output/figures/smoothness_bankloan1.eps", replace

}


if `figure_08' {

*------------------------------------------------------------------------------
* 1. SETUP AND DATA PREPARATION
*------------------------------------------------------------------------------

use AllIndiaII, clear

* Clean data
keep if abs(stdlog_i) < 4                    // Remove outliers beyond 3 SD
drop if log_i == .                           // Remove missing interest rates
drop if B_Ratio_Dist_Ex == .                 // Remove missing running variable

* Create household ID based on loan date for inflation adjustment
drop _merge
merge m:1 b14_q3 using CPI_index, keepusing(CPI_index)
drop if _merge != 3
drop if B_Ratio_Dist_Ex == .

* Inflation-adjust loan amount
gen loan_amount = b14_q5 / CPI_index * 100

* Treatment indicator (above threshold)
gen above = (B_Ratio_Dist_Ex > 0)

*------------------------------------------------------------------------------
* 2. GLOBAL BANDWIDTH CALCULATION
*------------------------------------------------------------------------------

* Estimate optimal bandwidth for main sample (2010-2013)
rdrobust log_i B_Ratio_Dist_Ex if inlist(b14_q3, 2010, 2011, 2012, 2013) & inlist(b14_q6,13,14), ///
    p(1) fuzzy(treatment) weights(Weight_SS) vce(cluster NSS_70_Codes) ///
    bwselect(mserd)

global h = e(h_l)  // Store bandwidth for subsequent analyses
di "Optimal bandwidth: $h"

drop _merge

*------------------------------------------------------------------------------
* 3. HETEROGENEITY BY WEALTH QUARTILES
*------------------------------------------------------------------------------

merge m:1 HHID using HH_wealth
drop if _merge != 3

* Create wealth quartiles
xtile quartile = value, nq(4)

* Save dataset with wealth information
save AllIndiaII_with_wealth, replace

*------------------------------------------------------------------------------
* 3A. RDD BY WEALTH QUARTILE (INTEREST RATE)
*------------------------------------------------------------------------------

matrix results = J(4, 3, .)  // quartile × (coef, CI_lower, CI_upper)

forvalues q = 1/4 {
    rdrobust log_i B_Ratio_Dist_Ex ///
        if inlist(b14_q3, 2010, 2011, 2012, 2013) ///
        & inlist(b14_q6, 13, 14) ///
        & quartile == `q', ///
        p(1) fuzzy(treatment) ///
        weights(Weight_SS) ///
        vce(cluster NSS_70_Codes) ///
        h($h)
    
    matrix results[`q', 1] = e(tau_cl)
    matrix results[`q', 2] = e(tau_cl) - 1.96 * e(se_tau_cl)
    matrix results[`q', 3] = e(tau_cl) + 1.96 * e(se_tau_cl)
}

* Prepare data for plotting
clear
matrix colnames results = coef ci_lower ci_upper
svmat results, names(col)

gen quartile = _n
label define quart_lbl 1 "1st" 2 "2nd" 3 "3rd" 4 "4th", replace
label values quartile quart_lbl

* WEALTH QUARTILE PLOT (INTEREST RATE)
twoway ///
    (rcap ci_upper ci_lower quartile, horizontal lcolor("217 85 130") lpattern(solid)) ///
    (scatter quartile coef, msymbol(O) mcolor("217 85 130")), ///
    xlabel(, angle(0) labsize(medium)) ///
    ylabel(1(1)4, valuelabel labsize(medium)) ///
    xline(0, lcolor(black) lpattern(dash)) ///
    xtitle("Estimated Coefficient", size(medium)) ///
    ytitle("Wealth Quartile Group", size(medium)) ///
    legend(off) ///
    plotregion(fcolor(white)) ///
    graphregion(fcolor(white) lwidth(0))

graph export "output/figures/rd_results_wealth_quartiles.pdf", replace
graph export "output/figures/rd_results_wealth_quartiles.eps", replace

*------------------------------------------------------------------------------
* 3B. HETEROGENEITY: HH HAS ANY BANK LOAN (BY WEALTH QUARTILE)
*------------------------------------------------------------------------------

use AllIndiaII_with_wealth, clear

* HH-level bank indicator: =1 if HH has at least one bank loan anywhere
bysort HHID: egen bank = max(inlist(b14_q6, 2, 3))

* Collapse to one row per HHID (keep variables needed for RD)
bysort HHID: keep if _n == 1

matrix results = J(4, 3, .)  // quartile × (coef, CI_lower, CI_upper)

forvalues q = 1/4 {
    rdrobust bank B_Ratio_Dist_Ex ///
        if inlist(b14_q3, 2010, 2011, 2012, 2013) ///
        & quartile == `q', ///
        p(1) fuzzy(treatment) ///
        weights(Weight_SS) ///
        vce(cluster NSS_70_Codes) ///
        h($h)

    matrix results[`q', 1] = e(tau_cl)
    matrix results[`q', 2] = e(tau_cl) - 1.96 * e(se_tau_cl)
    matrix results[`q', 3] = e(tau_cl) + 1.96 * e(se_tau_cl)
}

* Prepare data for plotting
clear
matrix colnames results = coef ci_lower ci_upper
svmat results, names(col)

gen quartile = _n
label define quart_lbl 1 "1st" 2 "2nd" 3 "3rd" 4 "4th", replace
label values quartile quart_lbl

* HH HAS ANY BANK LOAN (BY WEALTH QUARTILE) PLOT

twoway ///
    (rcap ci_upper ci_lower quartile, horizontal lcolor("217 85 130") lpattern(solid)) ///
    (scatter quartile coef, msymbol(O) mcolor("217 85 130")), ///
    xlabel(, angle(0) labsize(medium)) ///
    ylabel(1(1)4, valuelabel labsize(medium)) ///
    xline(0, lcolor(black) lpattern(dash)) ///
    xtitle("Estimated Coefficient", size(medium)) ///
    ytitle("Wealth Quartile Group", size(medium)) ///
    legend(off) ///
    plotregion(fcolor(white)) ///
    graphregion(fcolor(white) lwidth(0))

graph export "output/figures/rd_results_wealth_quartiles_bl_prob.pdf", replace
graph export "output/figures/rd_results_wealth_quartiles_bl_prob.eps", replace


*------------------------------------------------------------------------------
* 3C. HETEROGENEITY BY LOAN TYPE
*------------------------------------------------------------------------------

use AllIndiaII_with_wealth, clear

* Create loan type variable
gen loan_type = b14_q8
label values loan_type b14_q8

matrix results = J(4, 3, .)  // Matrix: type × (coef, CI_lower, CI_upper)

forvalues q = 1/4 {
    rdrobust log_i B_Ratio_Dist_Ex if inlist(b14_q3, 2010, 2011, 2012, 2013) & ///
        (loan_type == `q') & inlist(b14_q6, 13, 14), p(1) fuzzy(treatment) weights(Weight_SS) ///
        vce(cluster NSS_70_Codes) h($h)
    
    matrix results[`q', 1] = e(tau_cl)                     // Coefficient
    matrix results[`q', 2] = e(tau_cl) - e(se_tau_cl)*1.96 // Lower 95% CI
    matrix results[`q', 3] = e(tau_cl) + e(se_tau_cl)*1.96 // Upper 95% CI
}

* Save results for plotting
clear
matrix colnames results = coef ci_lower ci_upper
svmat results, names(col)

gen n = _n
label define type_lbl 1 "short pledged" 2 "short non-pledged" 3 "medium" 4 "long"
label values n type_lbl

* HETEROGENEITY BY LOAN TYPE PLOT

twoway (rcap ci_upper ci_lower n, horizontal lcolor("217 85 130") lpattern(solid)) ///
       (scatter n coef, msymbol(O) mcolor("217 85 130")), ///
       xlabel(, angle(0) labsize(medium)) ///
       ylabel(1(1)4, valuelabel angle(.1) labsize(medium)) ///
       xline(0, lcolor(black) lpattern(dash)) ///
       xtitle("Estimated Coefficient", size(medium)) ///
       ytitle("Loan type", size(medium)) ///
       legend(off) plotregion(fcolor(white)) graphregion(fcolor(white) lwidth(0))

graph export "output/figures/rd_results_by_type.pdf", replace
graph export "output/figures/rd_results_by_type.eps", replace

*------------------------------------------------------------------------------
* 3D. HETEROGENEITY -  LOAN TYPE & PROB(BEING A BANK LOAN)
*------------------------------------------------------------------------------

use AllIndiaII_with_wealth, clear

* Create loan type variable
gen loan_type = b14_q8
label values loan_type b14_q8

* HH-level bank indicator: =1 if HH has at least one bank loan anywhere
gen bank = inlist(b14_q6, 2, 3)

matrix results = J(4, 3, .)  // Matrix: type × (coef, CI_lower, CI_upper)

forvalues q = 1/4 {
    rdrobust bank B_Ratio_Dist_Ex if inlist(b14_q3, 2010, 2011, 2012, 2013) & ///
        (loan_type == `q'), p(1) fuzzy(treatment) weights(Weight_SS) ///
        vce(cluster NSS_70_Codes) h($h)
    
    matrix results[`q', 1] = e(tau_cl)                     // Coefficient
    matrix results[`q', 2] = e(tau_cl) - e(se_tau_cl)*1.96 // Lower 95% CI
    matrix results[`q', 3] = e(tau_cl) + e(se_tau_cl)*1.96 // Upper 95% CI
}

* Save results for plotting
clear
matrix colnames results = coef ci_lower ci_upper
svmat results, names(col)

gen n = _n
label define type_lbl 1 "short pledged" 2 "short non-pledged" 3 "medium" 4 "long"
label values n type_lbl

* LOAN TYPE & PROB(BANK LOAN TAKE-UP) PLOT

twoway (rcap ci_upper ci_lower n, horizontal lcolor("217 85 130") lpattern(solid)) ///
       (scatter n coef, msymbol(O) mcolor("217 85 130")), ///
       xlabel(, angle(0) labsize(medium)) ///
       ylabel(1(1)4, valuelabel angle(.1) labsize(medium)) ///
       xline(0, lcolor(black) lpattern(dash)) ///
       xtitle("Estimated Coefficient", size(medium)) ///
       ytitle("Loan type", size(medium)) ///
       legend(off) plotregion(fcolor(white)) graphregion(fcolor(white) lwidth(0))

graph export "output/figures/rd_results_by_type_bl_prob.pdf", replace
graph export "output/figures/rd_results_by_type_bl_prob.eps", replace


*------------------------------------------------------------------------------
* 3E. HETEROGENEITY BY YEAR
*------------------------------------------------------------------------------

use AllIndiaII_with_wealth, clear

matrix results = J(4, 3, .)  // Matrix: year × (coef, CI_lower, CI_upper)

forvalues q = 1/4 {
    local year = `q' + 2009  // 2010, 2011, 2012, 2013
    
    rdrobust log_i B_Ratio_Dist_Ex if inlist(b14_q3, `year') & inlist(b14_q6,13,14), ///
        p(1) fuzzy(treatment) weights(Weight_SS) ///
        vce(cluster NSS_70_Codes) h($h)
    
    matrix results[`q', 1] = e(tau_cl)                     // Coefficient
    matrix results[`q', 2] = e(tau_cl) - e(se_tau_cl)*1.96 // Lower 95% CI
    matrix results[`q', 3] = e(tau_cl) + e(se_tau_cl)*1.96 // Upper 95% CI
}

* Save results for plotting
clear
matrix colnames results = coef ci_lower ci_upper
svmat results, names(col)

gen n = _n
label define year_lbl 1 "2010" 2 "2011" 3 "2012" 4 "2013"
label values n year_lbl

* HETEROGENEITY BY YEAR PLOT

twoway (rcap ci_upper ci_lower n, horizontal lcolor("217 85 130") lpattern(solid)) ///
       (scatter n coef, msymbol(O) mcolor("217 85 130")), ///
       xlabel(, angle(0) labsize(medium)) ///
       ylabel(1(1)4, angle(1) valuelabel labsize(medium)) ///
       xline(0, lcolor(black) lpattern(dash)) ///
       xtitle("Estimated Coefficient", size(medium)) ///
       ytitle("Year when loan taken", size(medium)) ///
       legend(off) plotregion(fcolor(white)) graphregion(fcolor(white) lwidth(0))

graph export "output/figures/rd_results_by_year.pdf", replace
graph export "output/figures/rd_results_by_year.eps", replace

*------------------------------------------------------------------------------
* 3F. HETEROGENEITY BY PURPOSE
*------------------------------------------------------------------------------

use AllIndiaII, clear

gen purpose = 1
replace purpose = 6 if inlist(b14_q11, 1, 2)
replace purpose = 5 if inlist(b14_q11, 3, 4)
replace purpose = 4 if inlist(b14_q11, 8)
replace purpose = 3 if inlist(b14_q11, 10)
replace purpose = 2 if inlist(b14_q11, 11)
tab purpose

matrix results = J(6,3,.) // Empty matrix to store results

// Loop through each quartile and extract estimates
forvalues q = 1/6 {
	rdrobust log_i B_Ratio_Dist_Ex if inlist(b14_q3, 2010,2011,2012,2013) & inlist(purpose, `q') & inlist(b14_q6, 13, 14), ///
        p(1) fuzzy(treatment) weights(Weight_SS) vce(cluster NSS_70_Codes) h($h)

    // Extract coefficient and CI
    matrix results[`q',1] = e(tau_cl)       // Estimated coefficient
    matrix results[`q',2] = e(tau_cl)-e(se_tau_cl)*1.96       // Lower 95% CI
    matrix results[`q',3] = e(tau_cl)+e(se_tau_cl)*1.96       // Upper 95% CI
}

// Convert matrix to dataset
clear
svmat results
rename (results1 results2 results3) (coef ci_lower ci_upper)

// Create a quartile variable
gen n = _n
label define n_lbl 6 "farm business" 5 "non-farm business" 4 "education" 3 "medical treatment" 2 "housing" 1 "other"
label values n n_lbl

twoway (rcap ci_upper ci_lower n, horizontal lcolor("217 85 130") lpattern(solid)) ///
       (scatter n coef, msymbol(O) mcolor("217 85 130")), ///
       xlabel(, angle(0) labsize(medium)) ///
       ylabel(1(1)6, angle(1) valuelabel labsize(medium)) ///
       xline(0, lcolor(black) lpattern(dash)) ///
       xtitle("Estimated Coefficient", size(medium)) ///
       ytitle("Purpose of loan", size(medium)) legend(off) ///
       plotregion(fcolor(white)) graphregion(fcolor(white) lwidth(0)) //note("N = 23,393")
	   
* Save the graph as a PDF
graph export "output/figures/rd_results_by_purpose.pdf", replace
graph export "output/figures/rd_results_by_purpose.eps", replace

*------------------------------------------------------------------------------
* 3G. HETEROGENEITY IN PROB(TAKING-UP BANK LOAN) BY PURPOSE
*------------------------------------------------------------------------------

use AllIndiaII, clear

gen purpose = 1
replace purpose = 6 if inlist(b14_q11, 1, 2)
replace purpose = 5 if inlist(b14_q11, 3, 4)
replace purpose = 4 if inlist(b14_q11, 8)
replace purpose = 3 if inlist(b14_q11, 10)
replace purpose = 2 if inlist(b14_q11, 11)
tab purpose

* bank indicator: =1 if it is a bank loan
gen bank = inlist(b14_q6, 2, 3)

matrix results = J(6,3,.) // Empty matrix to store results

// Loop through each quartile and extract estimates
forvalues q = 1/6 {
	rdrobust bank B_Ratio_Dist_Ex if inlist(b14_q3, 2010,2011,2012,2013) & inlist(purpose, `q'), ///
        p(1) fuzzy(treatment) weights(Weight_SS) vce(cluster NSS_70_Codes) h($h)

    // Extract coefficient and CI
    matrix results[`q',1] = e(tau_cl)       // Estimated coefficient
    matrix results[`q',2] = e(tau_cl)-e(se_tau_cl)*1.96       // Lower 95% CI
    matrix results[`q',3] = e(tau_cl)+e(se_tau_cl)*1.96       // Upper 95% CI
}

// Convert matrix to dataset
clear
svmat results
rename (results1 results2 results3) (coef ci_lower ci_upper)

// Create a quartile variable
gen n = _n
label define n_lbl 6 "farm business" 5 "non-farm business" 4 "education" 3 "medical treatment" 2 "housing" 1 "other"
label values n n_lbl

twoway (rcap ci_upper ci_lower n, horizontal lcolor("217 85 130") lpattern(solid)) ///
       (scatter n coef, msymbol(O) mcolor("217 85 130")), ///
       xlabel(, angle(0) labsize(medium)) ///
       ylabel(1(1)6, angle(1) valuelabel labsize(medium)) ///
       xline(0, lcolor(black) lpattern(dash)) ///
       xtitle("Estimated Coefficient", size(medium)) ///
       ytitle("Purpose of loan", size(medium)) legend(off) ///
       plotregion(fcolor(white)) graphregion(fcolor(white) lwidth(0)) //note("N = 23,393")
	   
* Save the graph as a PDF
graph export "output/figures/rd_results_by_purpose_bl_prob.pdf", replace
graph export "output/figures/rd_results_by_purpose_bl_prob.eps", replace

*------------------------------------------------------------------------------
* 3H. HETEROGENEITY BY LOAN AMOUNT
*------------------------------------------------------------------------------

use AllIndiaII, clear

drop if b14_q5 == .

* Merge CPI data for inflation adjustment
drop _merge
merge m:1 b14_q3 using CPI_index, keepusing(CPI_index)
drop if _merge == 2
drop _merge

* CPI-adjusted loan amount (levels)
gen loan_amount = b14_q5 / CPI_index * 100

* Create wealth quartiles
xtile quartile = loan_amount, nq(4)

matrix results = J(4, 3, .)  // quartile × (coef, CI_lower, CI_upper)

forvalues q = 1/4 {
    rdrobust log_i B_Ratio_Dist_Ex ///
        if inlist(b14_q3, 2010, 2011, 2012, 2013) ///
        & inlist(b14_q6, 13, 14) ///
        & quartile == `q', ///
        p(1) fuzzy(treatment) ///
        weights(Weight_SS) ///
        vce(cluster NSS_70_Codes) ///
        h($h)
    
    matrix results[`q', 1] = e(tau_cl)
    matrix results[`q', 2] = e(tau_cl) - 1.96 * e(se_tau_cl)
    matrix results[`q', 3] = e(tau_cl) + 1.96 * e(se_tau_cl)
}

* Prepare data for plotting
clear
matrix colnames results = coef ci_lower ci_upper
svmat results, names(col)

gen quartile = _n
label define quart_lbl 1 "1st" 2 "2nd" 3 "3rd" 4 "4th", replace
label values quartile quart_lbl

* LOAN AMOUNT PLOT (INTEREST RATE)
twoway ///
    (rcap ci_upper ci_lower quartile, horizontal lcolor("217 85 130") lpattern(solid)) ///
    (scatter quartile coef, msymbol(O) mcolor("217 85 130")), ///
    xlabel(, angle(0) labsize(medium)) ///
    ylabel(1(1)4, valuelabel labsize(medium)) ///
    xline(0, lcolor(black) lpattern(dash)) ///
    xtitle("Estimated Coefficient", size(medium)) ///
    ytitle("Loan Amount Group", size(medium)) ///
    legend(off) ///
    plotregion(fcolor(white)) ///
    graphregion(fcolor(white) lwidth(0))

graph export "output/figures/rd_results_loan_amount.pdf", replace
graph export "output/figures/rd_results_loan_amount.eps", replace

*------------------------------------------------------------------------------
* 3I. HETEROGENEITY -  LOAN AMOUNT & PROB(BEING A BANK LOAN)
*------------------------------------------------------------------------------

use AllIndiaII, clear

drop if b14_q5 == .

* Merge CPI data for inflation adjustment
drop _merge
merge m:1 b14_q3 using CPI_index, keepusing(CPI_index)
drop if _merge == 2
drop _merge

* CPI-adjusted loan amount (levels)
gen loan_amount = b14_q5 / CPI_index * 100

* Create wealth quartiles
xtile quartile = loan_amount, nq(4)

* bank indicator: =1 if it is a bank loan
gen bank = inlist(b14_q6, 2, 3)

matrix results = J(4, 3, .)  // quartile × (coef, CI_lower, CI_upper)

forvalues q = 1/4 {
    rdrobust bank B_Ratio_Dist_Ex ///
        if inlist(b14_q3, 2010, 2011, 2012, 2013) ///
        & quartile == `q', ///
        p(1) fuzzy(treatment) ///
        weights(Weight_SS) ///
        vce(cluster NSS_70_Codes) ///
        h($h)
    
    matrix results[`q', 1] = e(tau_cl)
    matrix results[`q', 2] = e(tau_cl) - 1.96 * e(se_tau_cl)
    matrix results[`q', 3] = e(tau_cl) + 1.96 * e(se_tau_cl)
}

* Prepare data for plotting
clear
matrix colnames results = coef ci_lower ci_upper
svmat results, names(col)

gen quartile = _n
label define quart_lbl 1 "1st" 2 "2nd" 3 "3rd" 4 "4th", replace
label values quartile quart_lbl

* LOAN AMOUNT PLOT PROB(BANK LOAN)
twoway ///
    (rcap ci_upper ci_lower quartile, horizontal lcolor("217 85 130") lpattern(solid)) ///
    (scatter quartile coef, msymbol(O) mcolor("217 85 130")), ///
    xlabel(, angle(0) labsize(medium)) ///
    ylabel(1(1)4, valuelabel labsize(medium)) ///
    xline(0, lcolor(black) lpattern(dash)) ///
    xtitle("Estimated Coefficient", size(medium)) ///
    ytitle("Loan Amount Group", size(medium)) ///
    legend(off) ///
    plotregion(fcolor(white)) ///
    graphregion(fcolor(white) lwidth(0))

graph export "output/figures/rd_results_loan_amount_bl_prob.pdf", replace
graph export "output/figures/rd_results_loan_amount_bl_prob.eps", replace

}